 % ================================         
        % simulating Aiyagari         
        % ================================
            disp('Simulating Aiyagari')
             
        tic
       csim0=nan(Nvill*(vss+1),T);
       ksim0=nan(Nvill*(vss+1),T);
       income01=nan(Nvill*(vss+1),T);
       kni=nan(Nvill*(vss+1),1);
       ki=nan(Nvill*(vss+1),1);
       ci=nan(Nvill*(vss+1),1);
       yi=nan(Nvill*(vss+1),1);

        for nn=1:Nvill
            for jj=1:vss+1
                        k=zeros(T,1);
        c=zeros(T,1);
        y=zeros(T,1);
        clag=zeros(T,1);
        Dc=zeros(T,1);
        mm=zeros(T,1);
        k(1)=kinit;
        mm(1)=Z(income0ind(1))+(1+r)*k(1);
        for j=2:T
            % this is capital choice last period - starting
            % capital today
              k(j)=max(interp1(m(:,income0ind((nn-1)*(vss+1)+jj,j-1)),kp,mm(j-1),'linear','extrap'),phi);
            % this is consumption last period
            clag(j)=(mm(j-1)-k(j));
            y(j)=Z(income0ind((nn-1)*(vss+1)+jj,j));
            c(j-1)=clag(j);
            mm(j)=Z(income0ind((nn-1)*(vss+1)+jj,j))+(1+r)*k(j);
        end
        kend=max(interp1(m(:,income0ind((nn-1)*(vss+1)+jj,j-1)),kp,mm(j-1),'linear','extrap'),phi);
        c(T)=(mm(j-1)-k(j));
        csim0((nn-1)*(vss+1)+jj,:)=c;
        ksim0((nn-1)*(vss+1)+jj,:)=k;
        income01((nn-1)*(vss+1)+jj,:)=y;
            end
        end
        
        % save meand of capital - this is total excess demand given ri
        kni(:)=mean(ksim0(:,Tinit:T),2);
        ki(:)=mean(ksim0(:,Tinit:T),2);
        ci(:)=mean(csim0(:,Tinit:T-1),2);
        yi(:)=mean(income01(:,Tinit:T),2);
        knfreq=zeros(Nstar,NAiy);
        for i=Tinit:T
            if   k(i)<=phi%kp(1)
                knfreq(income0ind(i),1)=knfreq(income0ind(i),1)+1;
            else
                for j=2:NAiy
                    if and(k(i)<=kp(j,1),k(i)>kp(j-1,1))
                        knfreq(income0ind(i),j)=knfreq(income0ind(i),j)+1;
                    end
                end
            end
        end
        Phifreqdis=knfreq(:,1)./sum(knfreq,2);
        Phifreq=sum(knfreq(:,1))./sum(sum(knfreq));
        
        knfreq1=zeros(1,100);
        kgrid=linspace(-5,25,99);
        
        for i=Tinit:T
            if   k(i)<=kgrid(1)
                knfreq1(1)=knfreq1(1)+1;
            else
                for j=2:99;
                    if and(k(i)<=kgrid(j),k(i)>kgrid(j-1))
                        knfreq1(j)=knfreq1(j)+1;
                    end
                end
                if   k(i)>kgrid(99)
                    knfreq1(100)=knfreq1(100)+1;
                end
            end
        end
        knfreq1=knfreq1./sum(knfreq1);
        KNfreq=knfreq1;
        toc
        
        % normalise frequencies and save
        knfreq=knfreq./sum(sum(knfreq));
        % do everytthing again if upper bound on the capital grid is binding
        if sum(knfreq(:,NAiy))>0.0005
            disp('Care, the upper bound on the capital grid is binding.')
            pause
        end
